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•  Thermal  management  in  PEM  fuel  cell  stacks  are  modeled  through  a  mathematical  model. 

•  Thermal  management  of  stacks  with  U  and  Z  configurations  are  compared  to  each  other. 

•  Increasing  coolant  flow  rate  and  manifold  sizes  will  increase  temperature  uniformity. 

•  The  Z  configuration  is  not  appropriate  for  stacks  with  very  small  manifold  sizes. 
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The  narrow  range  of  operating  temperature  and  small  temperature  differences  between  the  stack  and 
the  ambient  have  made  the  fuel  cell  thermal  management  as  one  of  the  key  factors  that  influence  the 
performance  and  durability  of  polymer  electrolyte  membrane  fuel  cell  (PEMFC)  stacks.  In  the  present 
study,  an  analytical  model  is  developed  to  investigate  coolant  flow  and  temperature  distributions  within 
a  PEMFC  stack.  The  model  consists  of  a  coolant  flow  distribution  submodel  and  a  thermal  submodel  for 
determination  of  coolant  mass  flow  distribution  between  different  cooling  flow  fields  of  the  stack  and 
the  temperature  distribution  within  the  stack,  respectively.  The  coolant  mass  flow  rate  and  the  tem¬ 
perature  distributions  in  stacks  with  U  and  Z  configurations  are  compared  with  each  other  using  the 
developed  model.  The  test  results  of  two  65-cells  stacks  are  presented  to  verify  the  simulation.  The 
results  indicate  that  the  Z  configuration  results  in  more  uniform  temperature  distribution  than  the  U 
configuration  in  low  values  of  the  manifold  cross  sectional  area.  However,  the  Z  configuration  cannot  be 
applied  in  the  stacks  with  very  small  manifold  sizes.  A  parametric  analysis  is  also  conducted  to  assess  the 
effects  of  some  parameters  on  the  temperature  distribution  in  a  stack. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

With  the  advantages  of  low  working  temperature,  high  power 
density,  quick  start  up,  transient  ability  and  low  emissions  features 
PEM  fuel  cells  is  considered  to  be  the  most  promising  candidate  for 
the  next  generation  power  source  for  transportation,  stationary, 
auxiliary  and  portable  applications  [1-4].  Despite  of  wide  re¬ 
searches  and  progressions  on  fuel  cells,  several  technical  barriers  to 
their  commercialization  still  exist,  especially  about  their  durability 


*  Corresponding  author.  Institute  of  Materials  and  Energy,  Iranian  Space  Research 
Center,  7th  kilometer  of  Imam,  Khomeini  Ave.,  P.O.  Box  81955-174,  Isfahan,  Iran. 
Tel.:  +98  31  33222428;  fax:  +98  31  33222446. 

E-mail  addresses:  amir.mec.iut@outlook.com  (A.  Amirfazli),  asghari@fuelcell.ir, 
asghari@me.iut.ac.ir  (S.  Asghari). 

http://dx.doi.org/10.1016/jjpowsour.2014.06.073 

0378-7753 /©  2014  Elsevier  B.V.  All  rights  reserved. 


and  cost.  There  is  a  significant  amount  of  heat  generation  in  the  fuel 
cell  stack  due  to  electrochemical  reactions  and  electrical  re¬ 
sistances.  The  heat  generated  is  comparable  to  electrical  power 
output  and  should  be  removed  efficiently  to  avoid  overheating  of 
the  components  and  guarantees  the  favorable  working  tempera¬ 
ture  range  for  the  current  PEMFCs  which  is  usually  in  the  range  of 
60-80  °C.  Improper  thermal  management  and  non-uniform  tem¬ 
perature  distribution  within  the  fuel  cell  stack  may  cause  electro¬ 
lyte  drying  (global  or  local)  or  electrode  flooding  which  both  lower 
the  fuel  cell  performance  [5-7  .  On  the  other  hand,  the  tempera¬ 
ture  difference  between  the  PEMFC  and  the  ambient  is  very  small  in 
comparison  with  the  internal  combustion  engines.  Thus,  a  proper 
thermal  management  for  PEMFC  stack  is  very  challenging,  partic¬ 
ularly  in  automotive  applications  which  require  stacks  with  high 
power  output  and  high  power  density. 
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Nomenclature 

fiem 

Reynolds  number  of  coolant  flow  in  a  manifold 

segment 

APst 

pressure  drop  due  to  splitting  or  combining  of  two 

fm 

friction  coefficient  in  a  manifold  segment 

streams  in  a  manifold  segment 

Uch 

coolant  flow  velocity  in  a  cooling  flow  field  channel 

A  Pc 

pressure  drop  due  to  divergence  or  convergence  of  two 

Qc 

coolant  flow  rate  of  a  cell 

streams  in  T  shape  section  of  a  cell 

Tlch 

cross  sectional  area  of  a  cooling  flow  field  channel 

A  Pf-h 

pressure  drop  due  to  friction  and  bending  in  a  cooling 

Z 

number  of  channels  in  cooling  flow  field  of  a  cell 

flow  field  channel 

APCell 

total  pressure  drop  in  cooling  flow  field  of  a  cell 

APfm 

pressure  drop  due  to  friction  in  manifolds 

Fr 

relaxation  coefficient 

fist 

hydraulic  resistance  related  to  splitting  or  combining 

Qcorr 

corrected  coolant  flow  rate  of  a  cell  before  relaxing 

of  coolant  flow  in  manifold  segments 

k^corr 

corrected  velocity  of  coolant  flow  in  a  cooling  flow  field 

Rc 

hydraulic  resistance  related  to  divergence  or 

channel 

convergence  of  coolant  stream  in  T  shape  section  of  a 

Q 

heat  generation  per  cell 

cell 

1/ 

cell  working  voltage 

APm 

total  pressure  drop  in  a  manifold  segment 

/ 

output  electrical  current  of  stack 

APen 

pressure  drop  due  to  contraction  of  coolant  flow  in 

Tm 

coolant  mean  temperature  of  inlet  and  outlet  section  of 

inlet  section  of  a  cell 

a  cell 

APex 

pressure  drop  due  to  expansion  of  coolant  flow  in 

<4 

surface  heat  flux  of  a  cooling  flow  field  channel 

outlet  section  of  a  cell 

Ts(x) 

surface  temperature  of  a  cooling  flow  field  channel 

fien 

hydraulic  resistance  related  to  contraction  of  flow  in 

IUT 

index  of  uniform  temperature 

inlet  section  of  a  cell 

Tlsa 

surface  area  of  a  cooling  flow  field  channel 

fiex 

hydraulic  resistance  related  to  expansion  of  flow  in 

Qc  h 

contribution  of  a  flow  field  channel  from  generated 

outlet  section  of  a  cell 

heat  in  its  related  cell 

N 

number  of  cells 

P 

perimeter  of  a  cooling  flow  field  channel 

M 

number  of  nodes  in  cooling  flow  field  channels 

cp 

specific  heat  capacity  of  coolant 

Um 

coolant  flow  velocity  in  a  manifold  segment 

X 

distance  from  the  channel  inlet  along  a  cooling  flow 

9 

coolant  density 

field  channel 

M 

coolant  dynamic  viscosity 

h 

convective  heat  transfer  coefficient 

L 

length  of  manifold  segment 

mch 

coolant  mass  flow  rate  in  a  cooling  flow  field  channel 

Am 

manifold  cross  sectional  area 

T 

average  temperature  of  a  stack 

Qm 

volumetric  flow  rate  of  coolant  in  a  manifold  segment 

Tch J 

temperature  at  inlet  section  of  a  cooling  flow  field 

fifm 

hydraulic  resistance  related  to  friction  in  a  manifold 

channel 

segment 

p- 

i 

o 

temperature  at  outlet  section  of  a  cooling  flow  field 

f^hm 

hydraulic  diameter  of  manifolds 

channel 

Essentially,  the  issue  of  thermal  management  in  fuel  cell  stack 
can  be  viewed  at  two  levels;  at  the  cell  level  to  ensure  proper 
membrane  hydration  and  thereby  ensure  good  conductivity  of  the 
membrane,  and  at  the  system  level  to  ensure  temperature  unifor¬ 
mity  along  the  stack  and  keep  the  stack  from  heating  up  [8  . 

A  number  of  studies  have  been  conducted  on  the  thermal 
management  issue  in  the  PEM  fuel  cells.  In  contrast  with  the 
extensive  single  cell  models  available  in  literature,  limited  number 
of  studies  is  available  for  the  modeling  and  simulation  of  cooling 
behavior  of  PEM  fuel  cell  stacks.  Most  of  the  studies  are  devoted  to 
modeling  of  cooling  channel  on  bipolar  plates,  investigation  of 
various  cooling  strategies  and  selection  of  an  appropriate  cooling 
method. 

Yu  and  Jung  proposed  a  thermal  management  strategy  of  a 
PEMFC  system  with  a  large  active  area  for  transportation  applica¬ 
tions  [9  .  Their  model  was  composed  of  different  sub-models  for 
the  water  transport  through  the  MEA,  the  electrochemical  reaction 
in  the  cathode  catalyst  layer  and  the  two-dimensional  temperature 
distribution  within  the  fuel  cell.  They  concluded  that  the  fuel  cell 
operating  temperature  can  be  effectively  controlled  by  adjusting 
cooling  fan  operation. 

In  another  study  by  Asghari  et  al.,  design  of  cooling  flow  field  as 
well  as  thermal  management  subsystem  of  a  5  kW  PEM  fuel  cell 
system  were  investigated  [7].  The  number  of  parallel  channels  in 
the  parallel  serpentine  flow  field  was  selected  as  the  design 
parameter  and  its  optimum  value  was  obtained  by  compromising 
between  the  minimum  pressure  drop  of  coolant  across  the  flow 


field  and  the  maximum  temperature  uniformity  within  the  bipolar 
plate. 

Lopez-Sabiron  et  al.  used  a  theoretical  ID  model,  based  on 
electrochemical  and  thermal  equations,  to  estimate  the  perfor¬ 
mance  of  open-cathode  PEM  fuel  cells  as  well  as  to  define  an  op¬ 
timum  cooling  system  and  channel  geometry  for  the  cathode  side 

[1]. 

Sasmito  et  al.  carried  out  a  computational  study  of  PEFC  stacks 
with  various  cooling  strategies.  It  was  shown  that  liquid  cooling 
yields  best  performance  among  the  alternative  designs  examined 

[10]. 

Park  and  Li  developed  a  non-isothermal  stack  model  to  analyze 
the  effects  of  flow  variance  and  temperature  distribution  on  the 
performance  of  a  PEMFC  stack  11  .  The  stack  model  consists  of  the 
flow  network  solver  for  pressure  and  mass  flow  distributions  for 
the  reactant  gas  streams  and  cooling  water,  and  the  heat  transfer 
solver  for  temperature  distribution  among  the  cells  in  the  stack.  It  is 
shown  that  the  effect  of  temperature  is  dominant  on  the  cell 
voltage  variance  when  the  flow  variance  is  small  for  sufficiently 
uniform  distribution  of  reactant  flow  among  the  cells  in  the  stack. 
Furthermore  they  concluded  that  flow  and  temperature  distribu¬ 
tion  have  a  different  influence,  and  a  judicial  matching  of  their 
distribution  can  provide  the  ideal  uniform  cell  voltage  distribution. 

As  a  rule  of  thumb,  PEM  fuel  cell  stacks  above  5  kW  should  be 
water  cooled,  those  below  2  kW  air  cooled,  with  a  decision  for 
stacks  in  between  being  a  matter  of  judgment  12,13  .  One  of  the 
problems  encountered  when  stacking  fuel  cells  is  a  non-uniform 
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coolant  flow  distribution  between  cooling  flow  fields  of  cells,  which 
causes  an  uneven  temperature  distribution  within  the  stack. 
Therefore,  geometrical  design  of  the  cooling  flow  field  and  mani¬ 
folds  are  important  issues  to  be  considered. 

In  the  present  paper,  an  analytical  model  including  two  sub¬ 
models  is  developed  to  study  the  coolant  flow  and  temperature 
distribution  within  a  water-cooled  PEM  fuel  cell  stack.  The  coolant 
flow  distribution  submodel  determines  the  distribution  of  coolant 
flow  between  different  cells  in  the  stack.  With  considering  heat 
generation  in  each  cell  and  using  the  results  of  the  flow  distribution 
submodel,  the  thermal  submodel  yields  the  temperature  distribu¬ 
tion  within  the  fuel  cell  stack.  The  U  and  Z  stack  configurations  are 
investigated  and  compared  to  each  other  through  the  analytical 
model.  The  results  of  the  thermal  submodel  are  verified  by  some 
experimental  data.  Finally,  a  parametric  study  is  conducted  to 
investigate  the  effect  of  some  main  parameters  on  the  temperature 
uniformity  in  the  fuel  cell  stack. 

2.  Analytical  model 

The  analytical  model  consists  of  two  main  submodels:  (1)  a 
coolant  flow  distribution  model,  which  is  used  to  determine  the 
mass  flow  distribution  of  coolant  between  different  cooling  flow 
fields  of  the  fuel  cell  stack;  and  (2)  a  thermal  model,  which  com¬ 
putes  the  coolant  temperature  increase  along  the  flow  field  and 
also  the  temperature  distribution  within  the  fuel  cell  stack. 


By  the  coolant  flow  distribution  submodel,  the  flow  distribution 
is  calculated  using  the  Bernoulli  equation  and  the  fuel  cell  stack  is 
modeled  as  a  network  of  hydraulic  resistances  in  parallel  and  in 
series;  describing  the  resistances  in  the  cooling  flow  fields  and  in 
the  manifolds  (Fig.  1).  As  indicated  in  Fig.  1,  two  flow  patterns  are 
considered  in  the  current  analysis:  U-configuration  and  Z- 
configuration. 

2.2.  Coolant  flow  distribution  submodel 

The  coolant  flow  and  pressure  distributions  are  determined 
using  an  iterative  approach.  A  systematic  algorithm  is  constructed 
based  on  the  assumption  that  the  flow  distribution  in  the  fuel  cell 
stack  can  be  modeled  by  means  of  hydraulic  resistances  in  the 
manifold  and  in  the  cooling  flow  fields.  In  other  words,  the  stack  is 
represented  as  a  network  of  hydraulic  resistances,  in  parallel  and  in 
series  as  shown  in  Fig.  1. 

In  the  following,  the  pressure  drops  included  in  the  model  are 
described  and  the  solution  algorithm  is  explained  in  detail.  The 
pressure  drops  included  in  the  analytical  model  are  divided  into 
four  classes: 

1.  The  pressure  drop  due  to  splitting  and  combining  of  two 
streams  in  the  inlet  and  outlet  manifolds  (APst-in,  APst-out) 

2.  The  pressure  drop  due  to  friction  of  flow  in  the  inlet  and  outlet 
manifolds  (APfm_in,  APfm 

-out) 
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Coolant  flow  inlet 


Coolant  flow  outlet 


I _ , _ I 

I 

Outlet  manifold 
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Coolant  flow  inlet 


Coolant  flow  outlet 


- 1 - 

Outlet  manifold 
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Fig.  1.  Network  of  hydraulic  resistances;  (a)  U  configuration  and  (b)  Z  configuration. 
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3.  The  losses  caused  by  deflection  of  the  main  stream  in  the 
manifold  into  or  out  of  the  cells  (APc_in,  APc_0Ut) 

4.  The  losses  related  to  the  cooling  flow  field  channels  which 
consist  of: 

a)  Pressure  drop  due  to  the  effects  of  friction  and  bending 

b)  Pressure  drop  caused  by  expansion  and  contraction  of  the 
main  stream  into  or  out  of  the  channels. 

2  A  A.  Pressure  drop  due  to  splitting  and  combining  of  two  streams 
in  the  inlet  and  outlet  manifolds 

The  shape  of  the  junctions  between  the  manifolds  and  the 
diffuser  part  in  the  cooling  flow  field  of  the  cells  are  approximated 
to  a  converging/diverging  T  shape  for  which  the  pressure  drop  re¬ 
lations  are  obtained  from  IDELCHICK’s  hydraulic  resistances 
handbook  [14]. 

These  pressure  drops  are  determined  by  using  the  following 
relations: 


^A^st-in  [*']  —  0.5  •  Pst-in  [*'] '  P '  Um~ in  [*'] 

Qm-in  M 


APfm_out  [*]  —  0.5-Pfm_out  [i\  •  p  •  Um-0ut  [l] 

Pfm-outM  =fm- out  x 


D 


hm-out 


(4a) 


(4b) 


J  64/Pem_out  [i]  <2100 

/m-out  =  |  o.316/j?em_out  [i]0'25  Rem  >  2100 
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2 A3.  The  losses  caused  by  deflection  of  the  main  stream  in  the 
manifolds  into  or  out  of  the  cells 

The  losses  due  to  the  deflection  of  main  stream  into  or  out  of  the 


Um- in  HI  = 


A 


m- in 


u, 


m- out 


[«]  = 


Qm-out[i\ 


A 


m- out 


Um-m\i}2 

(la) 

cells  are  correlated  to  the  flow  and  area  ratios.  This  pressure  loss 
can  be  determined  from  the  following  equations: 

(lb) 

APc-in  [i\  =  0.5  Rc_in  [i]  -  p -  L/ch  [i]2 

(5) 

P'  Qm-out  [i\ 2 

(2a) 

APc-out  [i]  =  0.5-J?c_out  [i]  •  P  •  Uch  [i]2 

(6) 

(2b) 

Uch  [i]  =  Qc .[I] 
zAh 

(7) 

where  Dm_in  [i]  and  Um.ou t  [i]  are  coolant  flow  velocities  in  ith 
segment  of  the  manifolds,  Pst-in  [i]  and  Pst-out  [i]  are  hydraulic  re¬ 
sistances  related  to  splitting  and  combining  of  coolant  flow  in  the 
inlet  and  the  outlet  manifolds,  Amjm  and  Am_out  are  the  cross 
sectional  areas  of  the  manifolds,  and  Qm-in  [i]  and  Qm-out  III  are  the 
volumetric  flow  rates  of  the  coolant  in  ith  segment  of  the  manifolds. 

The  resistance  coefficients  due  to  splitting  and  combining  of  the 
fluid  streams  (Pst-in,  Rst-out)  are  related  to  the  flow  ratios  between 
the  flow  in  the  manifolds  and  in  the  cells.  These  relations  are  shown 
in  able  1.  In  the  cases  where  explicit  relations  for  hydraulic  re¬ 
sistances  are  not  available,  a  correlation  is  obtained  using  curve 
fitting  of  data  available  in  IDELCHICK’s  hydraulic  resistances 
handbook  [14  . 


where  Qc  [i]  is  the  coolant  flow  rate  in  the  ith  cooling  flow  field  of 
the  stack,  Uc h  [i]  is  the  coolant  flow  velocity  in  each  channel  of  the 
cooling  flow  field,  z  is  the  number  of  channels  in  the  cooling  flow 
field,  Ach  is  the  cross  sectional  area  of  each  channel  and  Pc_in  [i]  and 
Rc- out  [i]  are  the  hydraulic  resistances.  As  shown  in  able  1,  these 
hydraulic  resistances  are  functions  of  the  area  ratio,  the  ratio  be¬ 
tween  cross  sectional  area  of  a  manifold  and  cross  sectional  area  of 
all  cooling  channels  in  each  cooling  flow  field,  and  the  flow  ratio. 

2.1.4.  Pressure  losses  related  to  the  cooling  flow  field 

As  shown  in  Eq.  (8a),  the  pressure  drop  related  to  the  cooling 
flow  field  is  mainly  caused  by  different  factors: 


2.1.2.  Friction  losses  in  the  manifolds 

The  coolant  flow  in  the  manifolds  is  simulated  as  a  flow  in  a 
conduit.  With  this  assumption,  the  pressure  drop  can  be  deter¬ 
mined  by  conventional  internal  fluid  flow  relations  [15].  The 
amount  of  pressure  drop  due  to  friction  in  ith  segment  of  the  inlet 
manifold  is  determined  based  on  the  following  correlations: 


APfm-in  [ 

i]  =  0.5  -  Rfm_in  [I]  •  P  •  Lrm_in  [i]2 

(3a) 

Pfm-in  I1’] 

=  /m-in[l]  x 

(3b) 

/m-inW  = 

r  64/ftem  in;/|  Rem  <  2100 

(3c) 

:  \  0.3 1  / Rem_jn  [i]0'25  Rem>  2100 

P^m-in  [*] 

_  pUm- in  [^Aim-in 

(3d) 

where  Pem_in  [i]  is  the  Reynolds  number,  fmjm  [i]  the  friction  coef¬ 
ficient  and  Dhm-in  and  p  are  hydraulic  diameter  and  dynamic  vis¬ 
cosity  of  the  coolant,  respectively. 

The  amount  of  pressure  drop  due  to  friction  in  ith  segment  of 
the  outlet  manifold  is  determined  using  the  following  correlations: 


APCeii  M  —  -APch  [i]  +  APc_jn  [i}\  +  APc_out  [i]  (8a) 

AP ch  [i]  =  A  Pf_b\i]  +  APen[i]  +  APex[l]  (8b) 


a)  Pressure  drop  due  to  friction  losses  in  the  flow  field  channels 
and  rotation  of  the  coolant  flow  through  bending  sections  of  the 
flow  field  channels  (A Pf-b) 

b)  Pressure  drop  due  to  contraction  and  expansion  of  the  coolant 
flow  through  T  shape  section  between  each  cell  and  the  corre¬ 
sponding  manifold  segment  (APen  +  APex) 

A  numerical  simulation  can  be  employed  to  characterize  the 
pressure  drop  due  to  friction  and  bending  losses  (A Pf-b)  in  the  cooling 
flow  field.  Therefore,  a  numerical  model  of  the  flow  field  is  developed 
and  run  for  a  range  of  coolant  flow  rates.  In  this  way  a  set  of  data 
including  different  coolant  flow  velocities  and  their  corresponding 
pressure  drops  of  the  coolant  across  the  flow  field  is  obtained.  The 
relationship  between  the  pressure  drop  and  the  flow  velocity  is 
determined  by  fitting  a  second  order  polynomial  to  these  data: 

A  Pf-b  [<]  =  AUch  [i]2  +  BUch  [i] 


(9) 


A.  Amirfazli  et  al.  /  Journal  of  Power  Sources  268  (2014)  533—545 


Table  1 

Hydraulic  resistance  correlations  of  convergence  and  divergence  of  coolant  flow  in  cells  and  manifolds  [14]. 
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where  A Pf-b  is  the  pressure  drop  of  the  coolant  across  the  flow  field 
of  ith  cell,  ych  is  the  mean  velocity  of  the  coolant  in  the  channels 
which  is  obtained  from  the  flow  rates,  and  A  and  B  are  constants 
which  are  determined  by  curve  fitting.  The  first  order  and  the 
second  order  terms  come  from  the  friction  losses  and  the  channel 
bending,  respectively. 

Finally,  the  pressure  drops  due  to  expansion  and  contraction  of 
the  coolant  flow  at  entrance  and  exit  section  of  each  cell  are 
determined  using  the  following  relations: 

APen  [i]  =  0.5  Ren-p  Uch  [i]2  (10) 


A Pex  [i]  =  0.5Rex-pUch  [l]2  (11) 

where  Ren  and  Rex  are  hydraulic  resistances  due  to  expansion  and 
contraction  of  the  coolant  flow  in  inlet  and  outlet  portion  of  each 
cooling  flow  field  channel  15]. 

2.1.5.  Solution  algorithm 

The  solution  algorithm  starts  with  an  initial  estimate  of  the 
coolant  flow  rate  through  the  flow  fields.  A  uniform  flow  distri¬ 
bution  between  different  flow  fields  is  assumed  for  the  first 
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iteration  step.  Then,  the  coolant  flow  rates  in  different  sections  of 
the  manifolds  are  determined.  The  local  velocities  in  the  manifolds 
and  in  the  channels  of  the  flow  fields  are  computed  knowing  the 
flow  rates.  The  different  pressure  drops  are  calculated  using 
equations  ( 1 ) — (11 ).  The  next  step  is  to  calculate  the  local  pressures 
in  the  inlet  and  outlet  manifolds.  By  assuming  that  the  pressure  at 
exit  port  of  the  stack  is  zero,  the  local  pressures  in  the  inlet  and 
outlet  manifolds  are  calculated.  The  pressure  drops  in  the  cooling 
flow  fields  can  now  be  corrected  by  knowing  the  local  pressures  in 
the  inlet  and  outlet  manifolds.  New  cell  flow  rates  are  then  calcu¬ 
lated  from  the  corrected  coolant  flow  velocities  in  flow  field 
channels  using  the  following  equation. 


[Qc[i]]j+i  =  [Qc\i]]j  +  ([QcorrM  -  QcWlj)^  (13) 

where  [Qc  [/]]/+ 1  is  the  new  under  relaxed  flow  rate  of  ith  cell,  j  is 
the  numerator  of  iteration  and  Fr  is  the  under  relaxation  coefficient. 

Once  the  new  flow  rates  are  obtained,  they  are  compared  with 
the  initial  guess  or  the  previously  calculated  cell  flow  rates.  The 
calculations  continue  with  the  new  obtained  set  of  cell  flow  rates 
until  the  differences  between  the  new  and  old  values  (i.e.  residuals) 
are  become  within  the  convergence  criteria.  Thus,  the  convergence 
criterion  of  coolant  flow  distribution  submodel  is  the  amount  of 
residuals  of  cell  coolant  flow  rates  between  two  consecutive  iter¬ 
ations  which  are  calculated  using  the  following  equation. 


Qcorr  [i]  —  h^corr  M-CzAch)  (12) 

where  Qcorr  U]  is  the  corrected  flow  rate  of  ith  cell,  UC0YY  [i]  is  the 
corrected  local  velocity  in  the  channels  of  the  flow  fields,  z  is  the 
number  of  channels  in  the  cooling  flow  field  and  Ac h  is  the  cross 
sectional  area  of  each  channel. 

Under  relaxation  of  the  cell  flow  rates,  is  also  applied  to  ensure 
stable  computation. 


i=N. 


Residual  =  £  [&[%-,  -  [Qcffllj 

1=1 


(14) 


The  final  output  of  the  coolant  flow  distribution  model  is  the 
coolant  flow  distribution  between  different  flow  fields  of  the  stack. 
A  MATLAB  code  was  developed  using  the  above  algorithm  to 
determine  the  coolant  flow  distribution  in  the  stack.  A  simplified 
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Fig.  2.  Simplified  solution  algorithm  of  the  coolant  flow  distribution  and  the  thermal  submodels. 
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scheme  of  the  solution  algorithm  for  coolant  flow  distribution 
submodel  can  be  seen  in  Fig.  2. 

2.2.  Thermal  submodel 

The  thermal  submodel  is  developed  to  determine  the  temper¬ 
ature  distribution  in  the  stack.  Within  the  framework  of  the  present 
study,  each  cooling  flow  field  can  be  considered  as  a  parallel  flow 
heat  exchanger  in  which  several  channels  of  rectangular  cross- 
section  extend  in  parallel  to  each  other  from  the  inlet  manifold  to 
the  outlet  manifold  and  the  heat  is  transferred  from  the  cathode 
side  of  membrane  electrode  assemblies  (MEA)  to  the  coolant 
streams  through  bipolar  plates.  The  analysis  is  based  on  the 
following  assumptions: 


its  related  cell  which  is  obtained  by  dividing  the  rate  of  heat  gen¬ 
eration  in  each  cell,  Q.  [i],  by  the  number  of  channels  in  the  flow 
field,  P  is  the  perimeter  of  cross-section  of  the  channel,  Asa  is  the 
surface  area  of  each  channel  which  is  equal  to  product  of  the 
parameter,  P,  and  the  channel  length,  TCh-i  is  the  coolant  tempera¬ 
ture  at  inlet  section  of  the  flow  field  channel  in  each  cell,  qch  [i]  is 
the  surface  heat  flux  which  is  obtained  by  dividing  the  parameter 
Qch  [*]  by  the  surface  area  of  the  channel,  x  is  the  distance  from  the 
channel  inlet  along  the  channel,  mch  is  coolant  mass  flow  rate  in  the 
channel  obtained  by  dividing  the  cell  mass  flow  rate  by  the  number 
of  channels  in  the  cooling  flow  field. 

The  specific  heat  capacity  of  the  coolant  is  assumed  to  be  a 
function  of  temperature  [18]. 


1.  Each  cell  in  the  stack  has  its  own  cooling  flow  field. 

2.  Heat  transfer  is  steady  state.  Thus,  heat  generated  in  each  cell  is 
removed  completely  by  the  coolant  and,  therefore,  the  tem¬ 
perature  distribution  in  the  cell  remains  unchanged  with  time. 

3.  The  heat  generated  by  the  electrochemical  reaction  in  the  cath¬ 
ode  side  of  each  cell  transfers  to  its  adjacent  cooling  flow  field 
stream  through  the  bipolar  plates  by  conductive  heat  transfer. 
The  reaction  rate  over  the  active  area  of  the  MEA  is  assumed  to  be 
uniform.  This  yields  an  approximately  uniform  heat  generation 
over  the  cell’s  active  area  and,  therefore,  there  is  a  constant  heat 
flux  at  the  walls  of  channels  in  the  cooling  flow  field.  By  this 
assumption,  the  temperature  distribution  can  be  non-uniform 
over  the  active  area  in  each  individual  cell  in  the  stack. 


The  rate  of  heat  generation  per  each  cell  in  the  stack  at  any 
working  voltage  can  be  determined  using  the  following  equation 

[12]: 

Q  M  =  (1-23  —  V[i])I  (15) 

where  Q  [i],V  [i]  and  /  are  the  rate  of  heat  generation  in  ith  cell, 
working  voltage  of  ith  cell,  and  output  electrical  current  of  the 
stack,  respectively. 


4.  Neglecting  the  entrance  effects,  the  coolant  flow  and  the  ther¬ 
mal  conditions  are  assumed  to  be  fully  developed. 

5.  Thermal  resistances  at  the  surface  of  the  cooling  channels  are 
negligible. 

6.  Viscous  dissipation  in  the  coolant  flow  is  negligible. 

7.  Only  heat  transfer  by  conduction  in  the  through-plane  direction 
from  the  MEA  to  the  cooling  flow  field  is  considered.  Net 
conductive  heat  transfer  in  the  in-plane  directions  in  the  bipolar 
plates  is  neglected. 

8.  The  free  convective  heat  transfer  from  the  surfaces  of  the  stack 
to  the  ambient  air  is  negligible  [7  . 

9.  The  amount  of  temperature  increment  at  each  manifold 
segment  is  calculated  under  an  iterative  procedure. 


Each  individual  channel  in  the  cooling  flow  field  is  modeled  as  a 
tube  of  rectangular  cross  section  extending  from  the  inlet  manifold 
to  the  outlet  manifold  in  parallel  with  other  channels  of  the  same 
flow  field.  The  temperature  distribution  in  the  bipolar  plate  along 
the  flow  channels  is  determined  using  governing  equations  of  in¬ 
ternal  flow  in  a  tube  with  rectangular  cross  section  and  constant 
surface  heat  flux  as  following  [16,17]: 
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where  h  is  the  convection  heat  transfer  coefficient,  qc h[i]  is  the 
contribution  of  each  flow  field  channel  from  the  generated  heat  in 
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where  Tm  [i]  is  mean  of  the  inlet  and  outlet  coolant  temperatures  in 
ith  cooling  flow  field: 
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The  inlet  coolant  temperature,  Tc h-/,  is  known  and  the  outlet 
coolant  temperature,  Tc h-0,  is  determined  by  the  following 
equation: 


Ph-O  [>] 


Qch  ffl 
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+  Ph-i  H 


(19) 


Since  the  specific  heat  capacity  of  the  coolant  varies  with  the 
mean  temperature  and  the  mean  temperature  is  a  function  of  cell 
outlet  temperature  which  itself  is  unknown,  thus  the  outlet  coolant 
temperature  and  the  specific  heat  capacity  in  each  cooling  flow 
field  should  be  determined  simultaneously  by  an  iterative 
procedure. 

A  concept  of  Index  of  Uniform  Temperature  (IUT)  was  proposed 
by  Chen  et  al.  to  evaluate  the  degree  of  uniform  temperature  profile 


Fig.  3.  PEM  fuel  cell  system  used  in  the  experimental  tests. 
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across  the  cooling  plates  19].  In  the  present  study,  a  modified  form 
of  this  index  is  introduced  as: 


IUT 


Em  Em  | m,j)  -  r 

N  xM 


Eiii  Em  T(i,j) 

NxM 


(20) 


where  N  is  the  number  of  cooling  flow  fields  in  the  stack,  M  the 
number  of  nodes  taken  along  the  cooling  channel  from  the  inlet 
manifold  to  the  outlet  manifold,  7  is  the  average  temperature  of  the 
stack  and  T  ( ij )  is  the  local  temperature  in  bipolar  plate  at  ith 
cooling  flow  field  and  jth  node.  With  regard  to  the  definition  of  this 
index  in  equation  (20),  its  value  decreases  as  the  degree  of  tem¬ 
perature  uniformity  within  the  stack  increases.  A  simplified  scheme 
of  the  solution  algorithm  for  the  thermal  submodel  can  be  seen  in 
Fig.  2. 


3.  Experimental 

In  the  current  research,  a  series  of  experimental  tests  were 
conducted  in  order  to  build  confidence  in  the  model.  The  tests  were 
performed  using  a  home-made  PEMFC  system  (Fig.  3).  The  PEMFC 
system  contains  2  stacks.  While  the  stacks  are  electrically  con¬ 
nected  in  series,  the  coolant  and  the  reactant  gases  are  delivered  to 
them  in  parallel.  One  of  the  stacks  has  a  U-configuration  flow 
pattern,  and  the  other  one  has  a  Z-configuration  flow  pattern.  Each 
stack  contains  65  cells  with  a  cell  active  area  equal  to  225  cm2.  The 
MEA  consists  of  NRE-211  membrane,  catalyst  layer  with  a  total  Pt 
loading  of  0.4  mg  cirT2  at  both  the  anode  and  cathode  sides  and 
SGL  carbon  cloth  with  micro  porous  layer  as  the  gas  diffusion  layer 
(GDL).  A  more  detailed  description  of  the  anode  and  the  cathode 
flow  fields  can  be  found  in  our  previous  work  [7  . 

The  cooling  flow  field  was  of  parallel  serpentine  pattern 
including  10  channels  with  cross  section  of  3  mm  x  2  mm.  The 
coolant  inlet  and  outlet  manifolds  were  of  the  same  shape  and  size 
with  cross  section  of  5  cm  x  3  cm. 

Fuel  and  oxidant  subsystems  deliver  hydrogen  and  air  at  a 
constant  stoichiometry  of  1.2  and  2.0,  respectively.  Hydrogen  and 
air  are  humidified  up  to  95%  relative  humidity  prior  to  entering  the 
fuel  cell  stacks.  Operating  pressures  of  the  anode  and  the  cathode 
sides  were  kept  constant  at  0.4  bar  (g)  and  0.3  bar  (g),  respectively. 
Unutilized  fuel  from  the  anode  exhaust  is  recirculated  back  to  the 
feed  stream  via  a  recirculation  pump. 

During  the  tests,  the  coolant  temperature  at  the  inlets  of  both 
stacks  was  kept  constant  in  the  49  °C,  with  the  total  coolant  flow 
rate  of  52  L  min-1.  Therefore,  each  stack  approximately  received 
26  L  min-1  of  the  total  coolant  flow  rate.  The  total  coolant  flow  rate 
and  inlet  pressure  of  the  stacks  were  measured  via  a  digital  flow 
meter  and  a  pressure  transducer,  respectively.  The  flow  meter  was 
gpi  TM075-N,  3/4  in  and  the  pressure  transducer  was  of  Wikai  type 
A- 10. 

A  series  of  temperature  measurement  tests  were  carried  out  on 
the  stacks.  Since  determination  of  temperature  distribution  within 
the  bipolar  plates  was  impossible  during  fuel  cell  stack  operation,  it 
was  decided  to  measure  the  temperature  distribution  along  the 
bipolar  plate  edges  instead.  Temperature  measurements  were 
made  using  a  “DT-9875  Thermal  Imager”  infrared  thermal  camera. 


4.  Results  and  discussion 

The  model  was  constructed  based  on  the  homemade  U  and  Z 
shape  stacks  to  analyze  the  thermal  behavior  of  these  stacks  and 
also  to  verify  the  model’s  predictions  by  the  experimental  results. 
In  other  words,  the  geometrical  characteristics  and  working  con¬ 
ditions  of  the  stacks  are  used  as  input  parameters  for  the  model. 


The  PEM  fuel  cell  stacks  considered  for  the  simulation  contain 
65  cells  with  225  cm2  active  cell  area.  The  inlet  and  the  outlet 
manifold  dimensions  are  assumed  to  be  3  cm  x  5  cm.  Unless 
otherwise  stated,  it  is  assumed  that  there  is  a  uniform  cell  voltage 
distribution  through  the  stack.  The  current  density  at  working 
condition  is  assumed  to  be  1.0  A  cm-2  at  0.6  V  per  each  cell.  The 
cooling  flow  field  is  of  modified  parallel  serpentine  shape.  Table  2 
lists  a  summary  of  input  parameters,  dimensions  and  properties 
of  the  fuel  cell  stacks  used  in  the  present  model. 

Another  input  parameter  of  the  model  is  the  flow  rate-  pressure 
drop  characteristic  of  the  cooling  flow  field  as  described  previously 
(Eq.  (11))  which  obtained  via  a  numerical  model.  The  numerical 
model  was  composed  of  10  channels  of  rectangular  cross-section  of 
2  mm  x  3  mm  extend  in  parallel  to  each  other  from  the  inlet 
manifold  to  the  outlet  manifold.  Inlet  and  outlet  boundary  condi¬ 
tions  of  the  model  were  assumed  to  be  the  mass  flow  inlet  and  the 
pressure  outlet,  respectively.  The  number  of  meshes  for  the  model 
was  totally  equal  to  1,446,025  hexahedral  cells.  A  series  of  nu¬ 
merical  simulations  were  performed  to  characterize  the  variations 
of  coolant  pressure  drop  across  the  cooling  flow  field  as  a  function 
of  mean  velocity  of  the  coolant  in  the  channels  (Eq.  (11 )).  The  values 
obtained  for  A  and  B  are  10,543  (Ns2  m-4)  and  1728.6  (Ns  m-3), 
respectively.  An  example  of  the  contour  plot  of  pressure  drops 
obtained  for  0.4  L  min-1  of  coolant  flow  rate  per  cell  is  depicted  in 
Fig.  4. 

Fig.  5  shows  the  overall  coolant  flow  distribution  between 
different  cells  in  the  U  and  Z-shape  stacks.  As  seen  from  Fig.  5,  there 
is  a  non-uniform  flow  distribution  in  both  stacks.  The  U-shape  stack 
shows  a  gradually  increasing  flow  distribution  along  the  stack 
while  the  flow  distribution  in  the  Z-shape  stack  is  of  parabola  type. 
The  range  of  variations  in  cell  flow  rates  in  the  Z-shape  and  the  U- 
shape  stacks  are  about  0.21  L  min-1  and  0.12  L  min-1,  respectively. 
The  standard  deviations  of  cell  flow  rates  are  computed  as  0.0619 
and  0.0325  for  the  U  and  Z  configurations,  respectively.  Therefore,  it 
seems  that  the  Z-shape  stack  has  a  more  uniform  coolant  flow 
distribution  than  the  U-shape  stack. 

Fig.  6  shows  contour  plots  of  temperature  distribution  in  the 
stacks  with  U  and  Z  configurations.  The  contour  plots  show  a 
symmetric  parabola  temperature  distribution  for  the  Z  configura¬ 
tion  and  an  increasing  temperature  distribution  for  the  U  config¬ 
uration.  The  temperature  is  minimum  near  the  stack  inlet  and 
increases  as  the  distance  from  the  inlet  increases.  There  is  a  uniform 
temperature  distribution  in  the  stacks  at  the  local  points  near  to  the 
inlet  but  the  temperature  non-uniformity  prevails  at  the  local 
points  far  away  from  the  inlet.  The  values  of  IUT  are  computed  as 


Table  2 

Parameters  and  properties  of  the  fuel  cell  stacks  used  in  the  present  simulations. 


Component 

Parameters 

Symbols  Value 

Cooling  flow  field 

Cooling  channel  dimensions 

^ch 

3  mm  x  2  mm 

Channel  length 

Asa  IP 

0.56  m 

Number  of  channels 

z 

10 

Number  of  turns 

3 

Constant  A  in  Eq.  9 

A 

10,543 

Constant  B  in  Eq.  9 

B 

1728.6 

Coolant 

Flow  rate  per  stack 

26  L  min-1 

Inlet  temperature 

49  °C 

Density 

P 

982  kg  ur3 

Dynamic  viscosity 

A 

464  x  10”3  N  s  ur2 

Conduction  coefficient 

655.3  N  m2  k1 

Stack 

Number  of  cells 

N 

65 

Cell  active  area 

225  cm2 

Manifold  dimensions 

Am 

5  cm  x  3  cm 

Manifold  length  for  each  cell 

L 

1  cm 

Cell  voltage 

V[i] 

0.6  V 

Stack  output  current 

I 

220  A 
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Fig.  4.  Contour  plot  of  pressure  distribution  in  the  cooling  flow  field. 


1.30  and  1.26  for  the  U  and  Z  configurations,  respectively.  Since  the 
value  of  IUT  for  the  Z  configuration  is  lower  than  that  of  U  config¬ 
uration,  the  Z  configuration  can  result  in  a  more  uniform  temper¬ 
ature  distribution  in  a  stack  with  the  above  mentioned  geometrical 
specifications. 

41.  Verification  of  the  simulation  results 

The  PEM  fuel  cell  system  was  run  under  the  same  conditions  as 
those  considered  for  the  simulation  (Table  2).  A  series  of  temper¬ 
ature  measurement  tests  were  carried  out  on  the  stacks.  As  stated 
before,  since  measurement  of  temperature  distribution  within  the 
bipolar  plates  was  not  possible  during  the  fuel  cell  stack  operation 


by  common  methods;  it  was  decided  instead  to  measure  the  tem¬ 
perature  distribution  along  the  vertical  and  horizontal  surfaces  of 
the  stacks. 

In  this  part  of  the  study,  the  non-uniformity  of  cell  voltages  was 
taken  into  account  in  the  simulations  in  order  to  increase  the 
precision  of  the  model.  Therefore,  experimentally  measured  volt¬ 
ages  of  the  cells  during  the  tests  were  used  as  the  input  parameters 
of  the  model.  In  Fig.  7  a  bar  graph  of  cells’  voltage  destribution  for 
the  U  and  Z  configurations  is  illustrated.  The  infrared  pictures  of  the 
two  stacks  under  operation  are  shown  in  Fig.  8.  The  inlet  manifolds 
are  at  the  bottom  and  the  exit  manifolds  are  at  the  top.  Further¬ 
more,  the  horizontal  direction  of  coolant  flow  in  the  stacks  is  shown 
by  the  arrows  in  Fig.  8.  As  seen  from  these  photos,  the  regions  near 
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Fig.  5.  Coolant  flow  distribution  in  the  stacks;  (a)  U  configuration  and  (b)  Z 
configuration. 


Fig.  6.  Contour  plot  of  temperature  distribution  within  the  stacks;  (a)  U  configuration 
and  (b)  Z  configuration. 
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to  the  inlet  manifolds  have  lower  temperatures  than  the  regions 
near  to  the  exit  manifolds.  Additionally,  in  the  stack  with  the  Z 
configuration,  temperature  is  higher  at  the  middle  region  of  the 
stack.  This  is  due  to  the  fact  that  the  cells  located  in  this  region 
receive  lower  coolant  flow  rates.  But  in  the  stack  with  U  configu¬ 
ration,  the  temperature  is  lower  at  the  cells  located  near  to  the  inlet 
port  of  the  stack,  because  these  cells  receive  higher  flow  rates  of  the 
coolant.  Moreover,  in  these  pictures,  some  hot  points  are  observed 
on  the  faces  of  the  stacks  which  lead  to  some  fluctuations  in  tem¬ 
perature  readings.  The  cause  of  forming  such  hot  points  is 
explained  in  the  next  paragraphs. 

In  Fig.  9,  the  temperature  variations  along  the  centerline  of  the 
upper  face  of  the  stack  close  to  the  exit  coolant  manifold  obtained 
from  experimental  measurements  are  compared  to  those  obtained 
by  the  simulation.  A  reasonably  good  agreement  between  the 
simulation  results  and  the  experimental  data  can  be  observed.  As 
seen  from  Fig.  9,  temperature  profiles  show  an  increasing  distri¬ 
bution  along  the  stack  for  the  U  configuration  and  an  approxi¬ 
mately  parabola  distribution  for  the  Z  configuration. 

The  fluctuations  seen  in  the  model’s  predictions  come  from  an 
uneven  experimentally  measured  voltage  distribution  between 
different  cells  in  the  stack  which  result  in  an  uneven  heat  genera¬ 
tion.  But  the  origin  of  high  fluctuations  in  the  experimental  results 
is  different.  The  experimental  data  have  been  extracted  from 


infrared  pictures  by  a  software.  Since  bipolar  plates  have  been 
fabricated  by  machining  of  composite  blanks,  there  are  some  color 
variations  on  bipolar  faces  and  edges.  These  color  variations  lead  to 
different  emissions  and  different  temperature  readings.  Therefore, 
high  fluctuations  in  experimental  data  are  not  real  and  just  come 
from  color  variations  on  the  bipolar  surfaces. 

It  is  worthy  of  notice  that  Pei  et  al.  have  also  found  such  a 
parabolic  temperature  distribution  in  their  Z  shape  stack  [20  .  They 
measured  the  temperature  distribution  experimentally  and 
concluded  that  the  temperature  differences  in  the  stack  exhibited  a 
parabola  distribution  at  different  current  densities. 

Fig.  10  shows  the  comparison  between  the  numerically  and 
experimentally  obtained  temperature  variations  along  the  vertical 
edges  of  two  certain  cells  in  the  U  shape  stack  from  region  near  to 
the  inlet  manifold  to  region  near  to  the  outlet  manifold.  One  of 
these  cells  is  located  in  the  middle  of  the  stack  (Fig.  10a)  and  the 
other  one  is  located  close  to  one  end  of  the  stack  where  coolant 
inlet  and  outlet  ports  are  placed  in  (Fig.  10b).  As  seen  from  this 
figure,  there  is  a  reasonably  good  agreement  between  the  simula¬ 
tion  results  and  the  experimental  data.  As  expected,  the  tempera¬ 
ture  in  any  point  on  the  edge  which  is  near  to  the  exit  manifold  is 
higher  than  that  of  its  corresponding  point  near  to  the  inlet 
manifold.  This  is  due  to  the  fact  that  the  coolant  temperature  at  the 
exit  region  of  cooling  flow  field  is  higher  than  that  of  the  inlet 
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Fig.  7.  Cells’  voltage  distribution;  (a)  U  configuration  and  (b)  Z  configuration. 
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Fig.  9.  Comparison  of  temperature  variations  along  the  centerline  of  the  upper  face  of 
the  stack  close  to  the  exit  coolant  manifold  from  the  simulation  and  experiments  at  the 
nominal  working  conditions;  (a)  U  configuration  and  (b)  Z  configuration. 


Fig.  8.  Infrared  pictures  of  the  stacks;  (a)  U  configuration  and  (b)  Z  configuration. 

region;  therefore,  temperature  of  the  bipolar  plate  at  this  region 
should  be  higher  in  order  to  maintain  uniform  heat  flux. 

As  discussed  above,  comparison  of  the  simulation  results  with 
the  available  experimental  data  shows  that  the  proposed  model  can 
reasonably  predict  the  thermal  response  of  a  PEM  fuel  cell  stack. 
However,  more  experimental  measurements  and  data  are  required 
to  improve  the  accuracy  of  the  model.  For  example,  designing  an 
experimental  setup  to  measure  the  coolant  flow  distribution  be¬ 
tween  different  cells  in  a  stack  can  be  advantageous  for  improving 
the  precision  of  the  hydraulic  resistance  relations  and  therefore  the 
accuracy  of  the  coolant  flow  distribution  submodel. 

4.2.  Parametric  analysis 

In  the  previous  sections,  it  was  stated  that  maintaining  uniform 
coolant  flow  and  temperature  distributions  within  a  PEM  fuel  cell 
stack  are  desired.  These  distributions  are  affected  by  some  pa¬ 
rameters  like:  size  of  the  inlet  and  outlet  manifolds,  flow  distri¬ 
bution  pattern  in  cells  and  size  of  the  flow  channels,  coolant  flow 
rate,  etc.  For  example,  sufficient  temperature  distribution  unifor¬ 
mity  can  be  achieved  through  an  enlarged  manifold  and  reduced 
flow  channel  sizes.  However,  excessively  small  flow  channels  can 
lead  to  excessive  pumping  power  required  to  drive  the  coolant  flow, 
in  addition  to  other  problems  such  as  manufacturing  difficulty. 
Optimal  cooling  of  a  stack  can  be  attained  by  optimizing  the  size  of 
the  manifold  and  flow  channels  as  well  as  the  coolant  flow  rate. 


A  parametric  study  is  conducted  to  investigate  the  effects  of  the 
manifold  size  and  the  coolant  flow  rate  on  the  temperature  dis¬ 
tribution  in  a  PEM  fuel  cell  stack. 

4.2.1.  Manifold  size 

The  cross  sectional  area  of  the  inlet  and  outlet  manifolds  of  a 
stack  has  a  considerable  effect  on  the  coolant  flow  distribution 
between  different  cells  in  the  stack.  Hence,  the  degree  of  temper¬ 
ature  uniformity  within  the  stack  is  affected  by  the  manifolds  sizes. 
As  a  simplifying  assumption,  the  inlet  and  outlet  manifolds  are 
assumed  to  be  the  same  size.  Unless  otherwise  stated,  the  input 
parameters  for  the  simulations  are  the  same  as  those  presented  in 
Table  2.  A  uniform  cell  voltage  distribution  is  also  considered 
through  the  stacks. 

In  Fig.  11,  the  IUT  is  compared  between  the  U  and  Z  configura¬ 
tions  for  the  various  cross  sectional  area  of  the  manifolds.  With 
regard  to  the  definition  of  IUT,  its  value  decreases  as  the  degree  of 
temperature  uniformity  within  the  stack  increases.  The  Z  configu¬ 
ration  results  in  more  uniform  temperature  distribution  than  the  U 
configuration  in  low  values  of  the  cross  sectional  area.  However,  at 
the  larger  values,  both  configurations  result  in  the  same  tempera¬ 
ture  uniformity.  It  can  be  explained  by  the  fact  that  the  pressure 
losses  in  the  manifolds  increase  as  the  cross  sectional  area  de¬ 
creases.  On  the  other  hand,  coolant  distribution  between  cells  be¬ 
comes  more  non-uniformed  by  increasing  pressure  losses  in  the 
manifolds  while  keeping  the  pressure  losses  in  the  cells  constant. 

Much  smaller  manifold  sizes  were  not  considered  in  the  current 
analysis,  because  during  the  analysis  of  the  Z  configuration  stack,  it 
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Fig.  10.  Temperature  profiles  along  vertical  edges  of  two  certain  cells;  (a)  cell  located 
in  the  middle  of  the  stack  and  (b)  cell  located  close  to  the  inlet  of  the  stack. 


was  revealed  that  the  very  small  manifold  sizes  lead  to  backflow  in 
some  middle  cells  due  to  excessive  pressure  drop  in  the  inlet  and 
outlet  manifolds.  The  resulting  backflow  causes  hot  coolant  in  the 
outlet  manifold  to  stream  into  the  cell  cooling  flow  field  and  finally 
pour  into  the  inlet  manifold.  Hence,  overheating  of  those  cells 
which  experience  the  backflow  is  inevitable.  The  occurrence  of 
backflow  in  the  case  of  small  manifold  sizes  was  also  confirmed  in  a 
special  pipe  flow  network  analysis  software  through  a  simplified 
model  of  Z  configuration  stack. 


4.2.2.  Effect  of  total  coolant  flow  rate 

Fig.  12  depicts  the  effect  of  coolant  flow  rate  on  the  temperature 
uniformity  in  the  stacks.  As  expected,  the  value  of  IUT  decreases  by 
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Fig.  12.  Effect  of  coolant  flow  rate  on  the  temperature  uniformity. 


increasing  the  coolant  flow  rate;  however,  the  rate  of  decrease  is 
larger  at  the  lower  flow  rates  and  manifold  cross  sectional  areas. 
This  trend  is  less  significant  in  the  Z  configuration  than  the  U 
configuration.  Regarding  to  uniform  temperature  distribution 
criteria,  it  is  beneficial  to  increase  the  coolant  flow  rate  as  much  as 
possible,  whereas  higher  coolant  flow  rate  leads  to  higher  parasitic 
power  due  to  coolant  pumping.  Thus,  a  trade-off  between 
maximum  temperature  uniformity  and  minimum  parasitic  losses 
due  to  excessive  coolant  pumping  determines  optimized  amount  of 
total  coolant  flow  rate.  Moreover,  in  the  whole  of  coolant  flow 
range,  the  IUT  in  the  Z  configuration  is  lower  than  that  of  the  U 
configuration  with  the  same  manifold  cross  sectional  area.  Thus, 
the  Z  configuration  results  in  a  more  uniform  temperature.  It 
should  be  noted  that,  the  smaller  manifold  sizes  were  not  consid¬ 
ered  for  Z  configuration,  as  they  causes  backflow  in  some  middle 
cells. 


5.  Conclusion 

In  the  present  study,  a  non-isothermal  mathematical  model  was 
introduced  for  investigation  of  thermal  management  in  PEM  fuel 
cell  stacks.  The  model  was  formed  by  two  main  submodels,  a 
coolant  flow  distribution  submodel  and  a  thermal  submodel.  Two 
common  stack  configurations  were  compared  to  each  other  with 
assistance  of  the  model.  The  results  depicted  that  the  degree  of 
temperature  uniformity  is  more  in  the  Z  configuration  than  the  U 
configuration.  Simulation  results  were  verified  by  comparing 
experimentally  measured  temperature  profiles  along  the  vertical 
and  horizontal  edges  of  the  bipolar  plates  with  those  obtained  by 
the  analytical  model.  A  parametric  analysis  was  also  conducted  to 
determine  the  effects  of  various  factors  on  the  stack  temperature 
uniformity.  Among  miscellaneous  parameters,  the  manifold  size 
and  the  coolant  flow  rate  were  selected  to  investigate  their  effects 
on  the  temperature  uniformity  in  the  stacks.  A  modified  form  of 
standard  deviation  was  also  introduced  as  an  index  of  uniform 
temperature.  In  addition,  the  following  results  were  obtained: 

1  The  degree  of  temperature  uniformity  within  the  stack  increases 
as  the  cross  sectional  area  of  the  inlet  and  outlet  manifolds  in¬ 
crease.  The  Z  configuration  results  in  more  uniform  temperature 
distribution  than  the  U  configuration  in  low  values  of  the 
manifold  cross  sectional  area.  However,  at  the  larger  values, 
both  configurations  result  in  the  same  temperature  uniformity. 

2  The  value  of  IUT  decreases  as  the  coolant  flow  rate  increases.  The 
rate  of  reduction  in  IUT  is  larger  at  the  lower  flow  rates.  On  the 
other  hand,  higher  coolant  flow  rate  leads  to  higher  parasitic 
power  due  to  coolant  pumping.  Thus,  optimized  amount  of  total 
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coolant  flow  rate  is  determined  by  compromising  between 
maximum  temperature  uniformity  and  minimum  parasitic  los¬ 
ses  due  to  excessive  coolant  pumping. 
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